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Abstract We propose a Multiscale BDDC for a class of saddle-point prob- 
lems. The method solves for both flux and pressure variables. The fluxes are 
resolved in three-steps: the coarse solve is followed by subdomain solves, and 
last we look for a divcrgcncc-frcc flux correction and pressure variables us- 
ing conjugate gradients with a Multilevel BDDC preconditioner. Because the 
coarse solve in the first step has the same structure as the original problem, 
we can use this procedure recursively and solve (a hierarchy of) coarse prob- 
lems only approximately, utilizing the coarse problems known from the BDDC. 
The resulting algorithm thus first performs several upscaling steps, and then 
solves a hierarchy of problems that have the same structure but increase in 
size while sweeping down the levels, using the same components in the first 
and in the third step on each level, and also reusing the components from the 
higher levels. Because the coarsening can be quite aggressive, the number of 
outer iterations can be kept small and the additional computational cost is 
significantly reduced due to the reuse of the components. We also provide the 
condition number bound and numerical experiments confirming the theory. 

Keywords Iterative substructuring • balancing domain decomposition • 
BDDC • multilevel methods • multiscale methods • saddle-point problems 

Mathematics Subject Classification (2000) 65F08 • 65F10 • 65M55 • 
65N55 • 65Y05 



Supported in part by the National Science Foundation under grant DMS-0713876, and by 
the Grant Agency of the Czech Republic GA CR 106/08/0403. Support from DOE/ASCR 
is also gratefully acknowledged. 



B Sousedik 

Department of Aerospace and Mechanical Engineering, University of Southern California, 
Los Angeles, CA 90089-2531, USA. E-mail: sousedik@usc.edu 

Institute of Thermomechanics, Academy of Sciences of the Czech Republic, Dolejskova 
1402/5, 182 00 Prague 8, Czech Republic. 

Part of the work has been completed while the author was a Research Assistant Professor at 
the Department of Mathematical and Statistical Sciences, University of Colorado Denver. 



2 



Bedfich Sousedik 



1 Introduction 

The Balancing Domain Decomposition by Constraints (BDDC), proposed in- 
dependently by Cros [4], Dohrmann [6], and Fragakis and Papadrakakis [11], 
is along with the Finite Element Tearing and Interconnecting - Dual, Primal 
(FETI-DP) method by Farhat et al. [8,9] currently one of the most advanced 
and popular methods of iterative substructuring. These methods have been 
derived by modifications of the BDD method by Mandel [17], and of the FETI 
method by Farhat and Roux [10], respectively. The relations between these 
two families of methods have been studied extensively by many analysts in 
the substructuring field cf., e.g., [1,16,18,20], and also [25]. The methods have 
been also extended to multiple levels: one can find multilevel extensions of 
the BDDC in [21,26,27,32,33] and of FETI in [14]. From our current point of 
view arc important extensions to saddle-point problems, such as to the Stokes 
problem [15,23,24] and to the flow in porous media. One of the first domain 
decomposition methods for mixed finite element problems were proposed by 
Glowinski and Wheeler [12]. Their Method II has been preconditioned using 
BDD by Cowsar et al. [3] and using BDDC by Tu [31]. This approach is some- 
times regarded as hybrid because the method iterates on a system of dual 
variables (as Lagrange multipliers) enforcing the continuity of flux variables 
across the substructure interfaces. However, in order to allow for a multilevel 
extension, we would like to retain the original primal variables, and therefore 
we find the recent work of Tu [29,34] to be more relevant for our approach. 

The multiscale finite element methods introduced by Hou and Wu [13] 
have evolved in the recent years into an active research area. The methods 
gained popularity due to the ability to model problems which contain many 
spatial scales such as composite materials or porous media. Their main idea is 
to construct multiscale finite element basis functions that are adaptive to the 
local property of the differential operator and provide a way to handle systems 
with extremely large numbers of degrees of freedom especially due to highly 
heterogeneous media. These methods also provide an important and interest- 
ing alternative to the empirical upscaling methods cf., e.g., an overview by 
Cushman et al. [5] . From this perspective our method can be viewed as a way 
of numerical upscaling via the coarse basis functions known from the BDDC. 
The construction of basis functions is fully decoupled and thus amenable to 
massively parallel computers. However, because the BDDC method and its 
multilevel version as such is now only a part of larger algorithmic framework, 
we have decided to call this new method as a Multiscale BDDC. 

In this paper, we propose the Multiscale BDDC for a class of saddle-point 
problems. Our starting point is the algorithm of Ewing and Wang [7], see also 
Mathew [22]. The basic idea is to solve for flux variables in three-steps: first 
we perform a coarse solve which is followed by independent subdomain solves 
with zero boundary conditions in the second step. In the third step, we look 
for a flux correction and pressures. Due to the design of the algorithm, the flux 
correction is divergence-free, and we can use conjugate gradients (CG, resp. 
PCG) with a Multilevel BDDC preconditioncr that preserves all of the iterates 
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in the divergence-free subspace. Applications of two-, resp. three-level BDDC 
in the third step of this algorithm have been studied by Tu in [29,34]. Also, 
one has to make a careful decision in the design of the coarse solve for the first 
step. A straightforward idea is to use the same, but "coarse" finite element 
discretization and a natural (linear) interpolation between the two meshes 
as considered in [22,29]. Alternatively, the coarse solve has been obtained by 
an action of the BDDC preconditioner on a carefully chosen vector by Tu 
in [30,34] and she has also numerically observed a very similar performance 
of the two choices [30, Section 4.8]. Obviously, we will favor the second idea. 
Next, noting that the coarse solve in the first step has the same structure 
as the original problem, wc can use the algorithm recursively, and solve a 
hierarchy of coarse solves only approximately. The resulting algorithm thus 
first creates a hierarchy of problems with similar structure scaling-up through 
the levels. Then this hierarchy is solved, while sweeping down the levels in a 
loop of outer iterations, using the same components in the first and the third 
step on each level, and also reusing the components from all of the previous 
(higher) levels. Because the coarsening can be quite aggressive, the number 
of outer iterations can be kept small and the additional computational cost 
is significantly reduced due to the reusing of components. Unlike some of the 
previous works, we do not use the global partially assembled matrices neither 
the change of variables which seem to require more work in the present settings. 

It is important to note that for the solution of closely related Stokes 
problem, the algorithm is reduced to step three because the solution itself is 
divergence-free. We also remark that the present approach is limited by a spe- 
cial choice of finite elements. In particular, we will work with the lowest-order 
Raviart-Thomas (RTO) elements that have piecewise constant basis functions 
for pressure variables. This is not the case when, e.g., Taylor-Hood elements 
are used and the BDDC preconditioned operator is no longer invariant on the 
divergence- free subspace [24] . Finally, we note that our framework allows for ir- 
regular mesh decompositions, heterogeneous coefficients possibly utilizing the 
adaptive approach as in [19,26], and also allows for a relatively straightforward 
extension into 3D. However, such extensions will be studied elsewhere. 

The paper is organized as follows. In Section 2 we introduce the model 
problem, in Section 3 wc introduce its mixed finite element discretization and 
recall the original algorithm of Ewing and Wang. In Section 4 we derive the 
two-scale version of this algorithm using the BDDC components. In Section 5 
we formulate the Multiscale BDDC method. In Section 6 we derive the condi- 
tion number bound for the model problem, and finally in Section 7 we report on 
numerical experiments with a particular application to fiow in porous media. 

Throughout the paper we find it more convenient to work with abstract 
finite dimensional spaces and linear operators between them instead of the 
space M" and matrices. The results can be easily converted to the matrix 
language by choosing a finite element basis. For a symmetric positive definite 
bilinear form a, we will denote the energy norm by ||w||^ = a (u, u). 
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2 Model problem 

Let be a bounded polygonal domain in M", n = 2. We are interested in a 
solution of the following scalar, second-order, elliptic problem 

- V • kVp = f in f2, (1) 

where fc is a symmetric, uniformly positive definite matrix with bounded coeffi- 
cients, the right-hand side / € L? {Q), subject to sufficiently smooth boundary 
data on dQ = Fd U Fm- The variable p will be called pressure. 
Introducing the so-called flux variable 

u = -kVp, (2) 

we may rewrite (1) as a first-order system, generally known as Darcy's problem, 

k^^u + Vp = Q inQ, 
V • u = / in 

p = qn on Fm, 
u-n = gD on Fd, 

where n is the unit outward normal of ]7, and for the the boundary conditions 
it holds that gN E H^/^ (Fn), and go E H'^^^ (Fd). 
Let us define a space 

H(r2; div) = {v e ; V • V G L'^in); v ■ n = on dH D Fd} , 

and denoting by Hf2 the characteristic size of i7, the norm is defined as 

ll«llH(r2;div) = 11^11^2(^2) + l|V ■ V||^,(^) . 

The weak form of the Darcy's problem is 

/ k^^u-vdx— / pV-vdx= / (7Arv-nds, Vv e H(i7; div), (3) 
J^2 J n jPn 

V-uqdx = -( fqdx, yqeL^{Sl). (4) 
Jf2 Jn 

We will consider the case Fm = 0, so that the right-hand side term in (3) will 
vanish, which requires the compatibility condition 

fdx+ [ gDds = 0, (5) 

and the pressure p will be determined uniquely up to an additive constant. 
We refer to the monographs [2,28] for additional details. 
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3 Mixed finite elements and basic algorithm 

Let U be the lowest order Raviart-Thomas finite element space with a zero 
normal component on df2 and Q be a space of piecewise constants with a zero 
mean on Q. These two spaces, defined on the triangulation Th oi H where 
h denotes the mesh size, are finite dimensional subspaces of H(i7;div) and 
L'^{f2), respectively, and they satisfy a uniform inf-sup condition, see [2]. 

In the mixed variational formulation of the Darcy's problem, cf. eq. (3)- (4), 
we would like to find a pair {u,p) G {U, Q) such that 



a{u,v) + b{v,p) =0, Vu e ?7, (6) 

VqeQ. (7) 

Let us split the domain fi into non-overlapping subdomains i = 1, . . . , A^. 
Accordingly, let us split the solution spaces as 

U^Uo+ (©,= iC/0 + C/corr, (8) 

Q = ®toQ^■ (9) 



The spaces Uq, Qq are obtained by considering subdomains for a moment 
as macroelements, and the spaces t/^, Qi, for i = 1, . . . ^ N, are obtained by 
a restriction from the global solution spaces U, Q. The introduction of the 
auxiliary space ?7corr C U is motivated by an observation that in general 

U 7^ Uo + {(BtlU^) , (10) 

because the fluxes on subdomain interfaces might not be constant. Also, in or- 
der to write Uj = ®fLiUi we require zero fluxes over the subdomain interfaces. 
Finally, to determine the pressure p uniquely, we will consider the component 
Po G Qoj which is constant in each subdomain i?^, to have a zero average over 
the whole domain fl, and the components pi £ Qi to have zero averages over 
the subdomain Qi and identically equal to zero in other subdomains. We note 
that we will take an advantage of this splitting, in particular because for all 
uj £ Ui and qq G Qq, it holds, by the divergence theorem, that 

biui,qo) = - [ (V-U/) qo = 0. (11) 

The following algorithm is due to Ewing and Wang [7], cf. also Mathew [22]. 
Algorithm 1 Find the solution {u,p) e {U,Q) of the problem (6)-(7) as 

N 

M = UO + ^ + Ucorr, 
i=l 



in the following three steps: Compute 
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1. the coarse component (uo,Po) £ (C/oiQo) by solving 

a{uo,vo) +b{vQ,pQ) ^0, Vwo e C/o, (12) 
b {uq, go) = (/, Qo) , Vgo G Qo- (13) 

Note that because Qo £ Q> general 

b{uo,q) ^ {f,q) , VqeQ. 

2. the substructure components {ui,pi) € (Ui, Qi) for i ~ 1, . . . , N from 

a{ui,Vi) + b{v.i,pi) -a(uo,Wi) , Vwj £ t/^, 

^dd i/ie computed solutions as 

u* = Uq + ^ 



N 

Ui 
1 



.Due to the correction in the second step, and with respect to (9), we obtain 

b{u\q)^{f,q) yqeQ. (14) 

On the other hand, from (10), in general u* ^ u. Therefore, we also need 
3. the correction Ucorr S Ucorr C U . Considering 

substituting into (6)- (7) and using (14), compute {ucorr,p) from 
a {ucorr, v) +b {v , p) = —a {u*,v) , Vu G [/, 

b{Ucorr,q) =0, Vq G Q. 

Remark 1 Because the pressure components po, pi computed in the first two 
steps are tested only against proper subspaces of U , we simply disregard tlicm. 

The application of the BDDC preconditioner for the computation of itcon- 
for the two-, resp. three-level BDDC method has been studied by Tu [29,34]. 
However, comparing (12)-(13) with (6)-(7), we see that in fact we can use the 
same algorithm recursively, with multiple scales, to solve for both uq and Ucorr- 
But first, let us reformulate Algorithm 1 as the two-scale BDDC method. 
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4 Two-scale BDDC 

We begin by introducing the substructuring components. Let i7 be decomposed 
into nonoverlapping subdomains f2i, i = 1,...,A'^, also called substructures, 
forming a quasi-uniform triangulation of 12 with the characteristic subdomain 
size H. Each substructure is a union of the lowest order Raviart-Thomas 
(RTO) finite elements with a matching discretization across the substructure 
interfaces. Let Fi = df2i\df2 be the interface of the substructure ili and let 
r = U^^Fi. Let us denote by T the set of all faces between substructures, i.e., 
in the present context the set of all intersections Fij = FiHFj, i ^ j such that 
meas{Fij) > 0. Note that with respect to our discretization we define only 
faces, but no corners (nor edges in 3D) known from other types of substruc- 
turing. Let us also slightly generalize setting of coefficient k by allowing for 
coefficients ki which are symmetric and uniformly elliptic in each subdomain 
separately. 

The first step in substructuring is typically reduction of the problem to 
interfaces. Because only the flux variables are associated with the interfaces, 
we will introduce additional spaces and operators related only to the space U. 
For this reason we can for the moment closely follow our previous work [21]. 

Let Wi be the space of flux finite element functions on a substructure fii 
such that all of their degrees of freedom on dfii n df2 are zero, and let 

W ^Wix---x Wn. 

Now U C W can be viewed as the subspace of all function from W continuous 
across substructure interfaces. Define Uj C U as the subspace of functions 
that arc zero on the interface F, i.e., the space of "interior" functions. Denote 
by P the energy orthogonal projection from W onto Ui, 

P : w £ W I — 5- vj G Ui : a {vj, zj) = a {w, zj) , Vz/ G Uj. 

Functions from {I — P)W, i.e., from the nuUspacc of P, are called discrete 
harmonic; these functions are a-orthogonal to Uj and they are energy minimal 
with respect to increments in [//. Next, let W be the space of all discrete 
harmonic functions that are continuous across substructure interfaces, that is 

W ^{I -P)U, (15) 

and in particular, 

U = Ui®W, Ui ±a W. (16) 

The BDDC method is a two-level preconditioner characterized by the se- 
lection of certain coarse degrees of freedom. In the present setting, the value 
of a coarse degree of freedom on a face will be taken as the average of the 
fine scale degrees of freedom on the same face. Let W C W he the subspace 
of all functions such that the values of any coarse degrees of freedom have a 
common value and vanish on df2. We assume that 



a is positive definite on W. 



(17) 
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Define Wn C W a,s the subspacc of all functions such that their coarse 
degrees of freedom between adjacent substructures coincide, and such that 
their energy is minimal, and let us also define Wa C W as the subspacc of all 
function such that their coarse degrees of freedom vanish. Clearly, functions in 
Wn are uniquely determined by the values of their coarse degrees of freedom, 
and 

W = Wa® Wn, and Wa -La Wn- (18) 

Let be a projection from W onto U, defined by taking some weighted average 
of corresponding degrees of freedom on substructure interfaces. Finally, let us 
define as before, cf. eq. (9), 

Q = Qo®Qu 

where Qq consists of constant functions in each subdomain, such that 

/ Oo = and / Q/ = 0, i = 
Jn Jfii 

We remark that the restriction to Qi can be enforced by "bubble" constraints, 
and the basis for Qo can be obtained from the one of Q with only a little extra 
work. In the following, we will also need a subspacc of balanced functions 

Wb^{v^W -.b go) = Vgo e Qo} , 

for which we get the equivalence, cf. [28, Section 9.4], and also [22, eq. (15)], 

c|IHIa<lkllH(r2;div)<C|l«llL y{v,q^)^{wB,Q^). (19) 

Next, observe that it is only required for u* to satisfy (14). In particular, we 
do not need the substructures to form the same discretization as on the finite 
element level. Instead, we can conveniently retain the algebraic framework of 
the BDDC method introduced above and use its coarse problem in place of 
the coarse solve in Step 1. Specifically, let us set Uq = Wn- We are now ready 
to take the second look at Algorithm 1 and formulate its first modification. 

Algorithm 2 (Two-scale BDDC) Find the solution {u,p) e {U,Q) of the 
problem (6)- (7) by computing: 

1. the coarse component eW solving for wq G Wn the system 

a{wo,vn)+b{vn,Po) = 0, VSn<^Wn, (20) 
b{wo,qo) = {f,qo), Vqo e Qo, (21) 

and applying the projection 



uo = Ewq. 
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2. the substructure components uj G Ui solving 

a{ui,vi) + b{vi,pi) = -a (uo,w/) , Vuj G Ui, 

b {ui, qi) = (/, qi) - b {uq, qj) , Vqj G Qi, 

and combining the solutions u* = uq + uj. 

3. the correction and the pressure {ucomP) G {U,Q) from 

a (Ucorr, v) +b {v,p) = — fl (u* , v) , Vu G U, 

b{ucorr,q) ^0, yqeQ. 

Specifically, use the PCG method with the two-level BDDC preconditioner 
defined in Algorithm 3. 

Finally, combine the three solutions as 

U = Uq + Uj + Ucorr- 

Note that we again disregard the pressures po £^nd pi from Steps 1 and 2 
as in Algorithm 1. The algorithm of the two- level BDDC preconditioner used 
in Step 3 is closely related to the original version for elliptic problems, cf. [21, 
Algorithm 11]. For completeness its version for saddle-point problems follows. 

Algorithm 3 (Two-level BDDC preconditioner) Define the preconditioner 

(r,0) G (f/',Q') I — > iu,p) G {U,Q) as follows: 

Compute the interior pre- correction (uj,pi) G {Ui,Qi) from 

a{ui,zi) -\-b{zi,pi) = {r,zi) , Vz/ G Uj, 
b{ui,qi)^0, yqieQi. 

Set up the updated residual 

rB&U', {rB,v) = {rB,v) ~[a{ui,v) + b{v,pi)], \/v £ U. 

Compute the substructure correction wa £ Wa from 

a{wA,ZA) +b{zA,PiA) = {rB,EzA) , Vz/i G Wa, 
b{ui,qi)=0, yqi^Qi. 

Compute the coarse correction {wn,Po) G (Wn,Qi3j from 

a{wn,zn) -\-b{zn,Po) = {rB,Ezn) , Vz/j G Wn, 
b{wn,qn) = 0, Vqo e Qo- 

Add the averaged corrections 

ub ^ E {wa + Wn) ■ 
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Compute the interior post- correction {vi,qi) G {Ui,Qi) from 



a{vi,zi) +h{zi,qi) 
b{vi,qi) 



a{uB,zi) , 
b{uB,qi) , 



Vzi £ Ui, 
Vqj G Qi. 



Apply the combined corrections 



u ~ uj + ub - vi, 
p^pi +Po- qi- 



Remark 2 The solve in the space Wa gives rise to independent problems on 
substructures and the global coarse problem in the space Wn is exactly the 
same as the one used in Step 1 of Algorithm 2. 

Clearly, one could implement Step 3 of Algorithm 2 by eliminating first 
the interior unknowns in (Ui,Qi), i.e., performing the static condensation, 

iteratively solving the problem in the spaces (w, Qo^ , and recovering the 

interiors after the convergence. This would remove the interior pre-, and post- 
corrections from Algorithm 3, and the output of the preconditioner would not 
be required to be divergence-free but only balanced, cf. [28, Setion 9.4]; such 
approach might be also more appealing from the practical point of view. For 
a proof that given a sufficient number of constraints, the PCG method with 
the two-level BDDC preconditioner is invariant on the space of balanced, resp. 
divergence- free functions see [29, Lemma 2] or Lemma 1 in the next section. 

Due to the splitting of the spaces, in particular by relaxing only the ve- 
locity components in the action of the preconditioner, and with respect to 
the equivalence of norms (19), we can conveniently use the a— norm in the 
following estimate, cf., e.g., [28, Lemma 9.10] for more detailed discussion. 

Theorem 4 ([29, Theorem 1]) The condition number k of the two-level 
BDDC preconditioner from Algorithm 3 satisfies the bound 



Remark 3 In [29, Lemma 8], the suprcmum was taken over the space (/ — P) Wb 
of discrete harmonic balanced functioTL Nevertheless, the bound remains the 
same by considering the larger space Wb, cf. [21, Remark 16] for details. 

In Algorithm 2, the coarse problem used in Steps 1 and 3 is solved exactly, 
and therefore becomes a bottleneck in the case of many substructures. In the 
next section we will suggest its further modification by using it recursively for 
Step 1 on a multiple of different scales, leading to the Multiscale BDDC. 



K < u! — max 



sup 

w£Wb 



\\{I-P)Ew\\ 
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5 Multiscale BDDC 

We extend Algorithm 2 to multiple scales by using it recursively for Step 1, 
leading to a multilevel decomposition, and introducing thus a loop of outer 
iterations with the size given by the number of different decomposition scales. 

The substructuring components from Section 4 will be denoted by an ad- 
ditional superscript ^ , as f2l , i = 1, . . -N^, etc., and called level 1. The level 1 
coarse problem solved in (20)-(21) will be called the level 2 problem. It has 
the same finite element structure as the original problem (6)- (7) on level 1, so 
we put W}j = and Qq = Q^- Level 1 substructures are level 2 elements and 
level 1 coarse degrees of freedom are level 2 degrees of freedom. Repeating this 
process recursively, level £ — 1 substructures become level i elements, and the 
level £ substructures are agglomerates of level £ elements. Level £ substructures 
arc denoted by i7|, i = 1, . . . , N^, and they are assumed to form a conforming 
triangulation with a characteristic substructure size H^. For convenience, we 
denote by the original finite elements and put H'-* = h. The interface 
on level £ is defined as the union of all level £ boundary nodes, i.e., nodes 
shared by at least two level £ substructures, and we note that C F^~^. 
Level £ — 1 coarse degrees of freedom become level £ degrees of freedom. The 
shape functions on level £ are determined by minimization of energy with re- 
spect to level £ ~ 1 shape functions, subject to the value of exactly one level £ 
degree of freedom being one and others level £ degrees of freedom being zero. 
The minimization is done on each level £ element (level £ — 1 substructure) 
separately, so the values of level £ — 1 degrees of freedom are in general dis- 
continuous between level £ — 1 substructures, and only the values of level £ 
degrees of freedom between neighboring level £ elements coincide. 

The development of the spaces on level £ now parallels the finite element 
setting in Section 4. Denote = Wff^. Let W^ be the space of functions on 
the substructure f2f, such that all of their degrees of freedom on dSlf n dH are 
zero, and let 

W'^ = W(x---x W^t. 

Now C can be viewed as the subspace of all functions from that are 
continuous across the interface F^. Define Uf C as the subspace of functions 
that are zero on F^, i.e., the functions "interior" to the level £ substructures. 
Denote by the energy orthogonal projection from onto Uj. Denote 
by the energy orthogonal projection from onto Uj. Functions from 
(/ — P^) W^, i.e., from the nuUspace of P^, are called discrete harmonic on 
level £; these functions are a-orthogonal to J7| and energy minimal with respect 
to increments in Uj. Next, let be the space of discrete harmonic functions 
that are continuous across substructure interfaces on level £, that is 

W' = {I- P') U', 

and in particular, 

= t/f ® W^, uj ±a W^- (23) 



12 



Bedfich Sousedik 



Let C be the subspace of all functions such that the values of any 
coarse degrees of freedom on level £ have a common valuejor all relevant level £ 
substructures and vanish on df2f D dfi. Define Wfj C as the subspace of 
all functions such that their level £ coarse degrees of freedom between adjacent 
substructures coincide, and such that their energy is minimal, i.e., 

W^n = {l -P^)'^h, (24) 

and for completeness let us also define C as the subspace of all func- 
tions such that their level I coarse degrees of freedom vanish. Clearly, functions 
in Wfj are uniquely determined by the values of their level £ coarse degrees of 
freedom, and 

Wi la W^n, W^^W^^® W^n- (25) 

Let be a projection from onto U^, defined by taking some weighted 
average on F^. Finally, let us define as before, cf. eq. (9), 

Q' = Q'o®Q\, (26) 

where Qg consists of constant functions in each level i substructure, such that 

/ Q^=0 and / Q^ = 0, ^ = l,...,iV^. 
Let us also first define, for levels ^ = 1,...,L — l,a hierarchy of balanced spaces 

= {wEW^ -.b {w, go) = Vqo e Qi} ■ 
We arc now ready to generalize the two-scale Algorithm 2 to multiple scales. 

Algorithm 5 (Multiscale BDDC) Find the solution {u^,p^) £ {U\Q^) 

of the problem (6)-(7) in the following steps: 

fore^l,...L-l, 

If £ = L — 1, set-up and solve (by a direct solve) the coarse problem for Uq~^ ■ 
Else, set-up level t -\- 1 decomposition and formulate the scale £4-1 problem. 

end 

for 1 = L-l,...l, 

Find the substructure components u\ G Uf 

and combine the two solutions u*'^ = Uq + u\. 
Find the correction u^orr- G and p^ € ( we keep p^ only if £ = 1), using 

the PCG with the Multilevel BDDC preconditioner from Algorithm 6. 
Combine the three solutions as = Uq + u\ + u^corr- 
If i > 1, set Uq~^ = (the coarse component for the next iteration) . 



end 
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We note that the first loop provides a natural approach of scaling-up 
through the levels. The Multilevel BDDC preconditioner used in Step 3 of 
Algorithm 5 consists of recursive application of the two-level BDDC precondi- 
tioner for the approximate solution of the hierarchy of the coarse problems that 
were pre-computed in Step 1. Even though the preconditioner differs only little 
from its original version for elliptic problems described in [21, Algorithm 17], 
we again include its saddle-point version here for completeness. 



Algorithm 6 (Multilevel BDDC preconditioner) Define the preconditioner 
(/,0) G (t/^',g^') ^ e (C/^QO as follows: 

for k ~ £,..., L — 1, 



r%eU'^\ {r%, v'^) ^ {r%,v'^)~ [a {u'},v'')+b{v\ Vv'^ G U\ 




rfe 
I ' 



(27) 
(28) 



Set up the updated residual 




rk 




rk 

77j 



(29) 
(30) 



If k~ L — 1, solve the coarse problem directly and set 



u 



L 



else, set = Wfl and set up the right-hand side for level fc + 1, 



end 



for A; = L — 1, . . . , 



Average the approximate corrections, 




(31) 
(32) 
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(33) 
(34) 



Apply the combined corrections, 



k 



A.* I A? 



u 



p 



p'i+Po-q'i- 



end 



In order to guarantee that the Multilevel BDDC preconditioner is invariant 
on the space of divergence-free functions, we will need the following: 

Assumption 7 Let there exist a sufficient number of constraints as coarse 
degrees of freedom, such that on any decomposition level £ it holds that 



Remark 4 In implementation, we can quite easily satisfy Assumption 7 by 
prescribing coarse degrees of freedom as averages over every face on every 
decomposition level. Then the values of coarse degrees of freedom of functions 
from the space are zero, and the values of coarse degrees of freedom for 
functions from the space Wfj for all (pairs of) adjacent substructures coincide. 

Lemma 1 Let Assumption 7 be satisfied. Then, for any w £ Wg , it also 
holds that (/ - P^) E^w £ and the output of the Multilevel BDDC pre- 
conditioner in Algorithm 6 is divergence-free. 

Proof First, let us show that given w £ W^, for all go £ Ql, i ^ 1, . . . , L — 1, 

biw,qo) = b{{l -P')E'w,qo) ^0. 

Using (25), w = w/i + wn- Then, by Assumption 7 and (24), we obtain 
b{{l~P') E'wA,qo)+b{{l-P') E'wn,qo)=biwn,qo)^0, 

which shows that (/ - P^) E^w & W^. Next, note (23) and (26). In particular, 
on any level £ and for any uj € Uf we have by (11) that 

b{wn +ui,qo) ^b{wn,qo) =0, Vqo e Qo' 

and using (34), the interior post-correction satisfies 

b{wn +ui,qi) = 0, yqieQi. 

Finally, using (27) and (29) implies that u'' is divergence-free. □ 



E^wn = Wn , Vtuyj G Wfj . 
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Thus with a careful choice of the initial guess, the output of the Multilevel 
BDDC preconditioner is divergence-free and by induction all the PCG iter- 
ates, which arc linear combinations of the initial error and the outputs of the 
preconditioner, stay in the divergence-free subspace. 

Finally, using the equivalence (19), Lemma 1, and also with respect to 
Remark 3, the following variant of [21, Lemma 20] carries over. 

Lemma 2 If for some uj^ > 1, 

||(/-pOii;V||' <c.^||u;^||;, Vw'eW'B, e = i,...,L-i, (35) 
then the Multilevel BDDC preconditioner satisfies k < Y[i=i'^^- 



6 Condition number bound for the model problem 

We will now derive a condition number bound for the model problem, using the 
lower bound derived by Tu [34] , which is limited to a geometric decomposition 
of the domain i7 on every decomposition level. In particular: 

Assumption 8 Each suhdomain Slf,£~Q,...,L — l and i = 1, . . . , is 
quadrilateral. The subdomains also form on every decomposition level i a quasi- 
uniform coarse mesh of the domain Q with a characteristic mesh size . 

Let us apply the abstract results from [21] to the particular case of the low- 
est order Raviart-Thomas dicrctization and our model problem. First, define 
a subspace Wfj^ of discrete harmonic balanced coarse basis functions as 



W. 



nB 



Let II 

^Ila(f2?) energy norm of a balanced function w G (. — 

1, . . . , i — 1, restricted to subdomain i7f , i — 1, . . . N^, and let ||w||„ be the 
norm obtained by piecewise integration over each i?^. To apply Lemma 2 
to our model problem, we need to generalize the polylogarithmic estimate 
from Theorem 4 to coarse levels. To this end, let /^"'"^ : —5- C/^+^ be an 

interpolation from the level i coarse degrees of freedom (i.e., level £+1 degrees 
of freedom) to functions in another space U^'^^ and assume that, for all levels 
£ = 1, . . . , L — 1, and level £ subdomains Slf, i = 1, . . . , N^, the interpolation 
satisfies for all w S W^g and for all f^^'^^ the equivalence 

with Cj/cf < const boimdcd independently of . . . , H^^^ . 

Remark 5 Since = /, the two norms are the same on ~ — . 
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For the three-level BDDC for saddle-point problems with the lowest-order 
Raviart-Thomas finite element discretization in two dimensions ,^he result of 
Tu [34. Lemma 5.5], can be written in our settings for all w G W}jg and for 
all f2f as 

P^^linf) ^ Mliof) < c\ \\l''w\\l(n-f) ' (37) 

where is an interpolation from the coarse degrees of freedom given by the 
averages over substructure faces, and C2/cJ < const independently of H/h. 

We note that the level 2 substructures are called subregions in [34] and 
that = I. The assumption (36), along with the equivalence of seminorms 
on a factor space provided an equivalence of norms [21, Lemma 22], which 
implies by induction the following: 

Lemma 3 ([21, Lemma 23]) For alii ^ 0, . . . , L - 1, and i = l,...,N^, 

(38) 

with C2/C1 < , independently of H^,. . . , H^'^^ . 

Next, using Lemma 3, we generalize the polylogarithmic estimate from 
Theorem 4 to coarse levels. 

Lemma 4 For all suh structuring levels £ = 1, . . . , L — I, 

||(/-PVV||^<C^(l + log^y||«;^||^ Vw'eWi. (39) 

Theorem 9 The Multilevel BDDC for the model saddle-point problem in 2D 
with face coarse functions satisfies the condition number estimate 

Proof The proof follows from Lemmas 2 and 4, with ui^ = ^1 + log j . 
□ 

Remark 6 For L = 3 we recover the estimate by Tu [34, Theorem 6.2]. 

Corollary 1 In the case of uniform coarsening, i.e. with / H^^^ = H/h 
and the same geometry of decomposition on all levels €=1,...L — 1, we get 



K < C^-^l + log H/hf^^'^^ 



(40) 
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7 Numerical experiments 

Numerical examples are presented for a Darcy's problem on a square domain 
in 2D discretized by the lowest order quadrilateral Raviart-Thomas finite ele- 
ments (RTO). A square domain was uniformly divided into substructures with 
fixed H^/H^~^ ratio on each level £. The boundary conditions did not allow 
any flux across the boundary. The right-hand side was given by a unit source 
and sink in two distant corners of the domain, so that the compatibility condi- 
tion (5) was satisfied. The method has been implemented in Matlab and for the 
preconditioned gradients we have used zero initial guess and stopping criterion 
for a relative residual tolerance of 10~^. The results for different coarsening 
ratios H^/H^~^ (the relative subdomain size) and varying number of outer it- 
erations given by the number of levels L, are reported in Table 1. For each L, 
there were L — 1 outer iterations £, i.e., £ = 1, . . . , L ~1, consisting of the three 
steps described in Algorithms 2 and 5. In the third step the flux correction 
was computed by PCG with the {£ + l)-level BDDC preconditioner. 

From the results in Table 1 we can observe that with increasing number of 
levels, the condition number grows as predicted by Theorem 9 and in particular 
by formula (40). Also, it appears that for a fixed number of levels the condition 
number grows only mildly with increasing H^/H^~^ ratio. 

Finally, we remark that numerical experiments with the three-level method 
carried out by Tu [34] indicated that the method is independent of jumps in 
coefficients, if they are aligned with the decomposition into subdomains on 
the top level. Because we feel that this prevents a practical use of the pro- 
posed method for a realistic reservoir simulations with randomly distributed 
coefficients with large jumps, we will address this issue in a separate study. 
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